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Abstract 

A set of weakly interacting spin-^ Fermions, confined by a harmonic oscillator po- 
tential, and interacting with each other via a contact potential, is a model system 
which closely represents the physics of a dilute gas of two-component Fermionic 
atoms confined in a magneto-optic trap. In the present work, our aim is to present a 
Fortran 90 computer program which, using a basis set expansion technique, solves the 
Hartree-Fock (HF) equations for spin-^ Fermions confined by a three-dimensional 
harmonic oscillator potential, and interacting with each other via pair-wise delta- 
function potentials. Additionally, the program can also account for those anharmonic 
potentials which can be expressed as a polynomial in the position operators x, y, and 
z. Both the restricted-HF (RHF), and the unrestricted-HF (UHF) equations can be 
solved for a given number of Fermions, with either repulsive or attractive interac- 
tions among them. The option of UHF solutions for such systems also allows us to 
study possible magnetic properties of the physics of two-component confined atomic 
Fermi gases, with imbalanced populations. Using our code we also demonstrate that 
such a system exhibits shell structure, and follows Hund's rule. 
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Catalogue Identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, 
N. Ireland 

Distribution format: tar.gz 

Computers : PC's/Linux, Sun Ultra 10/Solaris, HP Alpha/Tru64, IBM/AIX 
Programming language used: mostly Fortran 90 

Number of bytes in distributed program, including test data, etc.: size of the 
gzipped tar file 371074 bytes 
Card punching code: ASCII 

Nature of physical problem: The simplest description of a spin ~ trapped sys- 
tem at the mean field level is given by the Hartree-Fock method. This program 
presents an efficient approach of solving these equations. Additionally, this 
program can solve for time-independent Gross-Pitaevskii and Hartree-Fock 
equations for bosonic atoms confined in a harmonic trap. Thus the combined 
program can handle mean-field equations for both the fermi and the bose par- 
ticles. 

Method of Solution: The solutions of the Hartree-Fock equation corresponding 
to the fermi systems in atomic traps are expanded as linear combinations of 
simple-harmonic oscillator eigenfunctions. Thus, the Hartree-Fock equations 
which comprises of a set of nonlinear integro-differential equation, is trans- 
formed into a matrix eigenvalue problem. Thereby, its solutions are obtained 
in a self-consistent manner, using methods of computational linear algebra. 
Unusual features of the program: None 



1 Introduction 



Over the last several years, there has been an enormous amount of interest 
in the physics of dilute Fermi gases confined in magneto-optic traps [T|2|3f4f5] . 
With the possibility of tuning the atomic scattering lengths from the repulsive 
regime to an attractive one using the Feshbach resonance technique, there has 
been considerable experimental activity in looking for phenomenon such as 
superfluidity, and other phase transitions in these systems [Tp] . This has led 
to equally vigorous theoretical activity starting from the studies of so-called 
BEC-BCS crossover physics[3], search for shell-structure in these systems[4], 
to the study of more complex phases [5]. As far as the spin of the fermions 
is concerned, most attention has been given to the cases of two-component 
gases which can be mapped to a system of spin-| atoms[3;4]. Therefore, in our 
opinion, a quantum-mechanical study of spin-| fermions moving in a harmonic 
oscillator potential, and interacting via a pair-wise delta function potential, 
can help us achieve insights into the physics of dilute gases of trapped fermionic 
atoms. 
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With the aforesaid aims in mind, the purpose of this paper is to describe a For- 
tran 90 computer program developed by us which can solve the Hartree-Fock 
equations for spin-| fermions moving in a three-dimensional (3D) harmonic 
oscillator potential, and interacting via delta-function potential. A basis set 
approach has been utilized in the program, in which the single-particle orbitals 
are expanded as a linear combination of the 3D simple harmonic oscillator ba- 
sis functions, expressed in terms of Cartesian coordinates. The program can 
solve both the restricted-Hartree-Fock (RHF), and the unrestricted Hartree- 
Fock (UHF) equations, the latter being useful for fermi gases with imbalanced 
populations. We would like to clarify, that as far as the applications of this 
approach to dilute Fermi gases is concerned, at present it is not possible to 
reach the thermodynamic limit of very large N, where N is the total number 
of atoms in the trap. However, we believe that by solving the HF equations for 
a few tens of atoms, one may be able to achieve insights into the microscopic 
aspects such as the nature of pairing in such systems. This program is an 
extension of an earlier program developed in our group, aimed at solving the 
time-independent Gross-Pitaevskii equation (GPE) for harmonically trapped 
Bose gases[6j. Thus the combined total program accompanying this paper can 
now solve for both Bose and Fermi systems, confined to move in a harmonic 
oscillator potential, with mutual interactions of the delta-function form. As 
with our earlier boson program, because of the use of a Cartesian harmonic 
oscillator basis set, the new program can handle trap geometries ranging from 
spherical to completely anisotropic, and it can also account for those trap 
anharmonicities which can be expressed as polynomials in the Cartesian coor- 
dinates. The nature of interparticle interactions, i.e., whether they are attrac- 
tive or repulsive, also imposes no restrictions on the program. We note that 
Yu et al. [4] have recently described a Hartree-Fock approach for dealing with 
two-component fermions confined in harmonic traps with spherical symmetry, 
employing a finite-difference-based numerical approach. However, we would 
like to emphasize that, as mentioned earlier, our approach is more general in 
that it is not restricted to any particular trap symmetry. Apart from describ- 
ing the program, we also present and discuss several of its applications. With 
the aim of exploring the shell-structure in trapped fermionic atoms, using 
our UHF approach we compute the addition energy for spherically trapped 
fermions for various particle numbers, and obtain results consistent with a 
shell-structure and Hund's rule. 

The remainder of the paper is organized as follows. In the next section we 
discuss the basic theoretical aspects of our approach. In section [3l we briefly 
describe the most important subroutines that comprise the new enlarged pro- 
gram. Section [4] contains a brief note on how to install the program and prepare 
the input files. In section [5] we discuss results of several example runs of our 
program for different geometries. In the same section, we also discuss issues 
related to the convergence of the procedure. Finally, in section [U we end this 
paper with a few concluding remarks. 
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2 Theory 



We consider a system of N identical spin-| particles of mass m, moving in 
a 3D potential with harmonic and anharmonic terms, interacting with each 
other via a pair-wise delta function potential. The Hamiltonian for such a 
system can be written as 

N N 

H = J2h(r t )+gJ25(v t -r J ), (1) 

i=l i>j 

where represents the position vector of i— th particle, g represents the 
strength of the delta-function interaction, and h(ri) denotes the one-particle 
terms of the Hamiltonian 

h(Ti) = + ~mfe 2 + tfv? + "IzD + V™ h ( Xi , y h (2) 

where u x , uj y and u z are the angular frequencies of the external harmonic po- 
tential in the x, y and z directions, respectively, and V anh (xi, y^ z,\) represents 
any anharmonicity in the potential. In order to parametrize the strength of 
the delta-function interactions, we use the formula g = 4nh2a in our program, 
where a is the s-wave scattering length for the atoms. Next we will obtain the 
RHF and the UHF equations for the system. 

Assuming that N = 2n, and that the many-particle wave function of the 
system can be represented by a single closed-shell Slater determinant, the 
RHF equations for the n doubly occupied orbitals {^(r), i = 1, . . . , n} of the 
system are obtained to be [7] 

(^ + ^El^(r)| 2 )^(r) = 6^(r). (3) 

3=1 

Similarly, for a system with n\ up-spin (a) fermions, and n 2 down-spin (/3) 
fermions (n\ + n 2 = N), the UHF equations for the up-spin orbitals can be 
written as|7j 

712 

(^ + ^El^f ) W| 2 )# ) (r) =e^ (Q) (r), (4) 
i=i 

where {ipi (r), i = l,...,ni} and {ipj(r), j = l,...,n 2 }, represent the 
occupied orbitals corresponding to the up and the down spins, respectively. 
Similar to the the UHF equations for the down-spin orbitals can be deduced 
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easily from Eq. (PJ. As in our earlier work on the bosonic systems[6], we 
adopt a basis-set approach and expand the HF orbitals in terms of the 3D 
Harmonic oscillator basis functions. This approach is fairly standard, and is 
well-known as the Hartree-Fock-Roothan procedure in the quantum chemistry 
community [7j. Thus, for the RHF case, the orbitals are expressed as 

^basis ^basis 

M r ) = C 3i®n xj ,n yj ,n zj {r) = G 3$n xi {x)<t>n yj {y)<t>n zj {z) , (5) 

3=1 3=1 



where Cji represents the coefficient corresponding to the j-th 3D harmonic 
oscillator basis function ^n X] ,n yj ,n zj {j), in the expansion of the i-th occupied 
orbital i^i(r), and Nbasis is the total number of basis functions used. Note 
that $n xj ,n vj ,n Z j( r ) i s itself a product of three linear harmonic oscillator eigen- 
functions of quantum numbers n X j, n y j, and n Z j. Therefore, a set of functions 
&n X j,ny j,n Z j( r )j f° r different values of n x j, n y j, and n Z j, will constitute an or- 
thonormal basis set, leading to an overlap matrix which is identity matrix. For 
the UHF case, the corresponding expansion for up-spin particles is 

4 a) (*) = Y^c^UuMKM^l (6) 

3=1 



from which the expansion for the down-spin particles can be easily deduced. 
Upon substituting Eqs. j5j) and ([6j), in Eqs. j3]) and (111), respectively, one 
can obtain the matrix forms of the RHF/UHF equations[7j. As outlined in 
our earlier work[6], numerical implementation of the approach is carried out 
in the so-called harmonic oscillator units, in which the unit of length is the 
quantity a x = J~^j- , and that of energy is fkJ x . The resulting matrix equation 
for the RHF case is 

FC(i) = £iC(i), (7) 



where represents the column vector containing expansion coefficients {Cji,j = 
1, . . . , N basis } ofipi, ii is the corresponding energy eigenvalue, and the elements 
of the Fock matrix F are given by 

F itj = EtSij + V™ h + g J i>m D k>l . (8) 

k,l=l 



Above 
Ei 



(9) 
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expressed in terms of aspect ratios 1 y = ^f- and 7* = ;jS V£j are the matrix 
elements of the anharmonic term in the confining potential, Dk,i = X^Li CkiCu 
is a density-matrix element, and Ji,j,k,i represents the 3D two-fermion repulsion 
matrix defined as 



Ji 



t ^Tl x j i TL x jTl x } c Tl x i ^ TlyiYlyjTlyfcYlyi '~^Tl z j i Tl z j7l z f i Tl z i . 



(10) 



Each one of the J matrices in Eq. f fTOl ). corresponding to the three Cartesian 
directions, can be written in the form 

oo 

/ d£<M£)<M£)M£)<M£), (11) 

— oo 



where £ is the corresponding Cartesian coordinate in the harmonic oscillator 
units. An analytical expression for J ni n 3 n k ni can be found in our earlier work[6]. 
In the UHF case, one obtains two matrix equations for the up/down-spin 
particles of the form 

fi(a)M°) = gWcg), (12) 



where the represents the Fock matrix for the up-spin particles given by 
Fij a) = Ei5 itj + Vjf 1 + gT,k}=x Ji,j,kM% 4 Q) is the energy eigenvalue, and 
D^i = YliliCfaC'ti, are the elements of the down-spin density matrix. We 
can easily deduce the form of the Fock equation for the down-spin particles 
from Eq. (|T2l) . In our program, HF Eqs. and f fT2l) are solved employing the 
self-consistent field (SCF) procedure, which requires the iterative diagonaliza- 
tion of the Fock equations[7]. 



3 Description of the program 



In this section we briefly describe the main program and various subrou- 
tines which constitute the entire module. As mentioned in the Introduction, 
the present program is an extension of our earlier program for bosons[6]. 
Thus the new program, which compiles as trap.x, can solve for: (a) time- 
independent Gross-Pitaevskii equation for bosons, and (b) Hartree-Fock equa- 
tions for fermions, confined in a trap. Therefore, most of the changes in the 
present program, as compared the earlier bosonic program, are related to its 
added fermionic HF capabilities. However, we have also tried to optimize the 
earlier bosonic module of the program wherever possible. A README file as- 
sociated with this program lists all its subroutines. Thus, in what follows, we 



6 



will describe only those subroutines which are either new (fermion related), 
or modified, as compared to the older bosonic code[6]. For an account of the 
older subroutines not described here, we refer the reader to our earlier work[6]. 
Additionally, with the aim of making the calculations faster, in the present 
code, we use the diagonalization routines of LAPACK library [8j, which re- 
quires the linking of our code to that library. Therefore, for this program to 
work, the user must have the LAPACK/BLAS program libraries installed on 
his/her computer system. The letter F or B has been included in parenthesis 
after the name of each subroutine to show whether the subroutine is useful for 
Fermionic or Bosonic calculations. If it is applicable for both, we denote this 
by writing BF. 



3.1 Main Program OSCL (BF) 



This is the main program of our package which reads the input data, dynam- 
ically allocates relevant arrays, and then calls other subroutines to perform 
tasks related to the remainder of the calculations. The main modification in 
this program, as compared to its earlier version[6], is that it now allows for 
input related to fermionic HF calculations. Thus, the user now has to specify 
whether the particles considered are bosons or fermions. If the particles con- 
sidered are fermions, one has to further specify whether the RHF or the UHF 
calculations are desired. For the case of UHF calculations, the user also needs 
to specify the number of up- and down-spin orbitals. Because of the dynamic 
array allocation throughout, no data as to the size of the arrays is needed from 
the user. The program will stop only if it exhausts all the available memory 
on the computer. There is one major departure in the storage philosophy in 
the present version of the code as compared to the previous one|6j in that now 
only the lower/upper triangles of most of the real-symmetric matrices (such 
as the Fock matrix) are stored in the linear arrays in the packed format. This 
not only reduces the memory requirements roughly by a factor of two, but 
also leads to faster execution of the code. 



3. 2 BECFERMI_ DRV (BF) 



This is the modified version of the old subroutine BEC_DRV, and is called 
from the main program OSCL. As its name suggests, it is the driver rou- 
tine for performing: (a) calculations of the bose condensate wave function for 
bosons, or (b) solving the RHF/UHF equations for fermions. Apart from allo- 
cating a few arrays, the main task of this routine is to call either: (a) routines 
BOSE_SCF or BOSE_STEEP depending upon whether the user wants to 
use the SCF or the steepest-descent approach meant for solving the GPE[6], 
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or (b) routines FERMI RHF or FERMI_UHF depending on whether the 
RHF or UHF calculations are to be performed. 

3.3 FERMI_RHF (F) 

This subroutine solves the RHF equations for the fermions in a trap using the 
SCF procedure, mentioned earlier. Its main tasks are as follows: 

(1) Allocate various arrays needed for the SCF calculations 

(2) Setup the starting orbitals. This is achieved by diagonalizing the one- 
particle part of the Hamiltonian. 

(3) Perform the SCF calculations. For this purpose, the two-particle integrals 
Ji,j,k,i (cf. Eq. (fTUl )) are calculated during each iteration[6j. If the user has 
opted for Fock matrix/orbital mixing, it is implemented using the formula 

= xmix R (i) + (1 - xmix) , 

where R^' is the quantity under consideration in the i-th iteration, and 
parameter xmix quantifying the mixing is user specified. Thus, if Fock 
matrix mixing has been opted, xmix specifies the fraction of the new 
Fock matrix in the total Fock matrix in the z-th iteration. If the user 
has opted for the orbital mixing, then each occupied orbital is mixed as 
per the formula above. The Fock matrix constructed in each iteration is 
diagonalized using the LAPACK routine DSPEVX|8], which can obtain a 
selected number of eigenvalues/eigenvectors of a real-symmetric matrix, 
as against traditional diagonalizers which calculate the entire spectrum 
of such matrices. We use DSPEVX during the SCF iterations to obtain 
only the occupied orbitals and their energies, thereby, leading to a much 
faster completion of the SCF process in comparison to using a diagonal- 
izer which computes all the eigenvalues/vectors of the Fock matrix. The 
occupied orbitals are identified according to the aufbau principle. 

(4) The total energy and the wave function obtained after every iteration are 
written in various data files so that the progress of the calculation can 
be monitored. This process continues until the required precision (user 
specified) in the total HF energy is obtained. 

3.4 FERMI_ UHF (F) 

In structure and philosophy this subroutine is similar to FERMI_RHF, ex- 
cept that its purpose is to solve the UHF equations for interacting spin-| 
fermions confined in a harmonic potential. Because there are two separate 
Fock equations corresponding to the up- and the down-spin fermions, the 
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computational effort associated with this subroutine is roughly twice that of 
routine FERMI RHF. 



3.5 BOSE_SCF (B) 

This subroutine aims at solving the time-independent GPE for bosons us- 
ing the iterative diagonalization approach, and was described in our earlier 
paper [6j. The diagonalizing routine which was being used for the purpose ob- 
tained all the eigenvalues and eigenvectors of the GPE, which is quite time 
consuming for calculations involving large basis sets. Since the condensate cor- 
responds to the lowest-energy solution of the GPE, using diagonalizing rou- 
tines which obtain all its eigenvalues and eigenvectors is wasteful. Therefore, in 
the new version of BOSE_SCF we now use the LAPACK|8] routine DSPEVX 
to obtain the lowest eigenvalue and the eigenvector of the Hamiltonian during 
the SCF cycles, leading to substantial improvements in speed. 



3.6 BOSE_ STEEP (B) 

This subroutine aims at solving the time-independent GPE for bosons using 
the steepest-descent method, and was also described in our earlier paper [6]. 
In this routine, the main computational step is multiplication of a trial vector 
by the matrix representation of the Hamiltonian. In the earlier version of the 
code, because the entire Hamiltonian was being stored in a two-dimensional 
array, we used the Fortran 90 intrinsic subroutine MATMUL for the purpose. 
However, now that we only store the upper triangle of the Hamiltonian in 
a linear array, it is fruitful to use an algorithm which utilizes this aspect. 
Therefore, we have replaced the call to MATMUL by a call to a routine called 
MATMUL_UT written by us. This has also lead to significant speed improve- 
ments. 



3.7 MATMUL _ UT (B) 

As mentioned in the previous section, the aim of this subroutine is to multiply 
a vector by a real-symmetric matrix, whose upper triangle is stored in a linear 
array. This routine is called from the subroutine BOSE_STEEP, and it utilizes 
a straightforward algorithm for achieving its goals by calling two BLAS|8j 
functions DDOT and DAXPY. 
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3.8 Plotting Subroutines (BF) 



We have also significantly improved the capabilities of the program as far as 
plotting of the orbitals and the associated densities is concerned. Now the or- 
bitals, or corresponding densities, can be computed both on one-dimensional 
and two-dimensional spatial grids, along user-specified directions, or planes. 
The driver subroutine for the purpose is called PLOT_DRV, which in turn 
calls the specific subroutines suited for the calculations. These subroutines 
are PLOT1D, and PLOT1DUHF for the one-dimensional plots, and 
PLOT 2D and PLOT2DUHF for the planar plots. The output of this 
module is written in a file called orb_plot .dat, which can be directly used in 
plotting programs such as gnuplot or xmgrace. 



4 Installation, input files, output files 

In our earlier paper, we had described in detail how to install, compile, and run 
our program on various computer systems[6]. Additionally, we had explained 
in a step-by-step manner how to prepare the input file meant for running the 
code, and also the contents of a typical output file[6|. Because, various aspects 
associated with the installation and running of the program remain unchanged, 
except for some minor details, we prefer not to repeat the same discussion. 
Instead, we refer the reader to the README file in connection with various 
details related to the installation and execution of the program. Additionally, 
the file 'input_prep.pdf' explains how to prepare a sample input file. Several 
sample input and output files corresponding to various example runs are also 
provided with the package. 



5 Calculations and Results 

In this section we report results of some of the calculations performed by our 
code on fermionic systems. We present both RHF and UHF calculations for 
various types of traps. Further, we discuss some relevant issues related to the 
convergence of the calculations. 

5.1 RHF Calculations: total energy convergence 

In this section our aim is to investigate the convergence properties of the total 
HF energy of our program with respect to: (a) number of particles in the trap, 
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Table 1 



Convergence of total HF energy (Ehf) for a spherically symmetric trap containing 
two particles, with respect to the size of the basis set, for various positive values 
(repulsive interactions) of the scattering length. Above, nmax is the maximum value 
of the quantum number of the SHO basis function in a given direction, and N^asis 
is the total number of basis functions corresponding to a given value of nmax. In 
some cases, Fock matrix mixing approach was used to achieve convergence. 

(b) symmetry of the confining potential, (c) nature and strength of interac- 
tions, and (d) number of basis functions employed in the calculations. As far 
as the number of particles is concerned, we have considered two closed-shell 
systems namely with two particles (N = 2), and with eight particles (N — 8). 
For N = 2 case, calculations have been performed for all possible trap geome- 
tries ranging from a spherical trap to a completely anisotropic trap. During 
these calculations, we have considered both attractive and repulsive interac- 
tions, corresponding to negative and positive scattering lengths, respectively. 
The magnitude of the scattering length (|a|) employed in these calculations 
ranges from 0.01a x to 0.8az. To put these numbers in perspective, we recall 
that in most of the atomic traps, a x « 1.0 /im, and for a two-component 6 Li 
trapped gas, the estimated value of the scattering length is anomalously large 
a ~ — 2160a |9j, where a is the Bohr radius. Thus, for this very strongly in- 
teracting system, the scattering length a ~ — 0.1 la x , is well within the range 
of the scattering lengths considered in these calculations. Therefore, the sys- 
tems considered here — ranging from weakly interacting ones to very strongly 
interacting ones — truly test our numerical methods. 

The results of our calculations are presented in tables CD— [5l For N = 2 sys- 
tem, we performed these calculations in order to understand the convergence 
behavior of the total energy with respect to the basis set size, with the goal 
of a high precision (six decimal digit convergence) in the total energy. Such 
high accuracy on larger systems will be computationally much more expensive, 
and, therefore, our aim behind the study of N = 8 system was to understand 
the role of number of particles on our results. The next larger closed-shell sys- 
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Table 2 



Convergence of total HF energy for a spherically symmetric trap containing two par- 
ticles, with respect to the size of the basis set, for various negative values (attractive 
interactions) of the scattering length o. Various symbols have the same meaning as 
in table [TJ 

tern will correspond to iV = 20, but we have not studied that here, because, 
in our opinion, such calculations will not lead to any newer insights into our 
approach. Next we discuss our results on these systems individually. 

With the aim of a more detailed exposition of the convergence behavior for 
repulsive and attractive interactions, for N = 2 system corresponding to an 
isotropic trap, we present our results for the positive and negative scattering 
lengths in separate tables CD and El For the rest of the cases, results for the 
attractive and the repulsive interactions are presented in the same tables. 
Upon examining our results for iV = 2 case (cf. tables CD— Q]), we conclude that 
for the case of repulsive interactions, calculations always exhibit convergence 
from above on E HF , with respect to the basis set size. In order to achieve 
six-digit accuracy for repulsive interactions, one needs to use relatively large 
basis sets, although a three-digit accuracy can be obtained using considerably 
smaller basis sets. However, quite expectedly, a drastically distinct convergence 
behavior is seen for the cases involving attractive interactions. It is obvious 
that for the attractive interactions, for sufficiently large scattering length, the 
HF method will not be applicable, and will exhibit instabilities because of pair 
formation. For relatively weaker attractive interactions, one again encounters 
convergence from above, as was the case for repulsive interactions. But, as the 
strength of the attractive interactions increases, the convergence with respect 
to the basis set size becomes more difficult to achieve, and for \a\ > 0.3a x 
(with a < 0), this property is completely lost, and the HF method begins to 
exhibit unstable behavior. 
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Table 3 



Convergence of total HF energy for a cylindrical potential (7^ = 1, j z = \/8) 
containing two particles, with respect to the size of the basis set, for various values 
of the scattering length a. Above, nxmax is the maximum value of the quantum 
number of the SHO basis function in x- and y— direction, nzmax is the same number 
corresponding to the ^-direction. Rest of the quantities have the same meaning as 
explained in the caption of table [TJ In some cases, Fock matrix mixing was employed 
to achieve convergence. 

Inspection of tables [3] and [4] reveals that for a given value of interaction 
length, the convergence requires the use of larger basis sets with increas- 
ing trap anisotropy, ranging from the perfectly spherical traps, to completely 
anisotropic traps. This behavior is expected for cases with aspect ratios 7^ and 
7 Z > 1, because the effective interaction constant in such cases g' = ^J^ffug > 

Upon examining our results for N = 8 case (c/. table [5|), we again see very 
monotonic convergence behavior for all calculations corresponding to repul- 
sive interactions, and note that the high accuracy in Ehf can be achieved 
with reasonably sized basis functions. However, as was the case for N = 2, 
completely different behavior is encountered when the interactions are attrac- 
tive. The calculations with a = — 0.05a x exhibit systematic convergence in 
E HF with the increasing basis set size, but for the case with a = — 0.1a x , no 
trend towards the convergence emerges, pointing again towards an unstable 



13 











a = — 0.3a x 


a = — 0.1a x . 


a = 0.1 a x 


a = 0.20a; 


a = 0Aa x 


nxmax 


nymax 


nZTflCLX 


^ * basis 


1 ' 1 1 1 


1 ' 1 1 1 






' '11 1 


9 
Z 


n 
u 


n 
U 


Q 
O 


4.00010O 


o.ooyoo4 


fi ^79^01 
O.o i ZOU1 


fi 71 9A7D 
0. 1 1Z4 / U 


I .oZoOOO 


9 
Z 


9 

z 


n 
u 


Q 

y 


/I ^9/19Q7 
4.0Z4Zy ( 


O.O/I 4ZU 


D.oDoy / D 


D.Doooyy 


1 .Zo0y40 


9 
z 


9 

z 


9 

z 


97 
Z 1 


A ^771 QzL 
4.0/1 iy4 


0.00 / 04 / 


R '} c ;897fi 
O.oOoZ / 


O.DOOoll 


/ . looyoo 


/I 


9 
Z 


9 
Z 


A ^ 


/i 9^7f;8n 

4.ZO / OoU 


O.ODZODo 


O.oO / OoO 


0.004oy0 


7 1 s^o^n 
/ . looyou 


/I 




9 
Z 


7^ 


A 1 /1 1 Q/lfi 
4. 14iy40 


c; c:cn979 
O.ODUZ / Z 


D.oooy / y 


0.00o04y 


7 18/1917 


4 


4 


4 


19^ 
1Z0 


4.U0Z01U 


O.OOo / 4y 


O.oOOOOO 


D.DOZOoU 


7 1 S99D1 
/ . loZZul 


D 


/| 


4 


1 7^; 

1(0 


q qi /|a«n 
O.y 14oOU 


0.00 1 4Uo 


O.oOOOUo 


D.oozooy 


7 1 89Dfi8 
/ . loZUOo 


u 


u 


u 


040 


^ 74Q' : i1 4 


O .OOVOiJ L 


U.OOU4UU 


u.uuzooy 


7 1 81 QD^ 


8 


8 


8 


819 


3.344656 


5.555972 


6.356391 


6.662373 


7.181859 


10 


10 


10 


1331 


2.746459 


5.555775 


6.356390 


6.662370 


7.181852 


12 


12 


12 


2197 


1.871496 


5.555707 


6.356390 


6.662369 


7.181851 



Table 4 



Convergence of total HF energy for an anisotropic potential = 2, 7 2 = 3) con- 
taining two particles, with respect to the size of the basis set, for various values of 
the scattering length. Above, nxmax, nymax, and nzmax represent the maximum 
values of the quantum number of the SHO basis function in x-, y— , and z— direc- 
tions, respectively. Rest of the quantities have the same meaning as explained in 
the caption of tabled! In some cases, Fock matrix mixing was employed to achieve 
convergence. 

behavior. 

Finally, in Fig. Q] we present the orbital density plots for the N = 2 case 
with both attractive and repulsive interactions, corresponding to a = ±0.2a x . 
The noteworthy point in the graph is the accumulation of the density at the 
center of the trap in case of attractive interactions, as compared to when 
the interactions are repulsive. With increasingly attractive interactions, this 
phenomenon becomes even more pronounced, possibly causing the instabilities 
in the HF approach. 



5.2 Unrestricted Hartree-Fock Calculations 

In this section we describe the results of our UHF calculations. If one performs 
a UHF calculation on a closed-shell system, one must get the same results as 
obtained by an RHF calculation. Similarly, the total energy and orbitals of 
a system with m up-spin and n down-spin particles should be the same as 
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Table 5 



Convergence of total HF energy (Ehf) for a spherical symmetric potential contain- 
ing eight particles, with respect to the size of the basis set, for various values of the 
scattering length. Different symbols above have the same meaning as explained in 
the caption of table [TJ In all the calculations presented above, Fock matrix mixing 
was used to achieve convergence. 




Figure 1. Density p(r) = 2\^i s (r)\ 2 plotted along the x-axis, obtained from RHF 
calculations on a two-particle system in an isotropic trap with a = 0.2a x (solid lines), 
and a = —0.2a x (dashed lines). Distance r is in harmonic oscillator units. 

that of a system with n up-spin and m down-spin particles. These properties 
of the UHF calculations can be used to check the correctness of the underly- 
ing algorithm. We verified these properties explicitly by: (a) performing UHF 
calculations on closed-shell systems with various scattering lengths and ge- 
ometries, and found that the results always agreed with the corresponding 
RHF calculations, and (b) by performing UHF calculations on various open- 
shell systems with interchanged spin configurations and found the results to 
be identical. Therefore, we are confident of the essential correctness of our 
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Figure 2. Calculated UHF values of addition energies (A / u(A r ) = fj,(N + 1) — fi(N)) 
of a spherical trap (in the units of hjj x ) with scattering length a = 0.01o x , plotted 
as a function of the particle number N, ranging from N = 1 to N = 21. 

UHF program, and in what follows, we describe its applications in calculat- 
ing the addition energy of fermionic atoms confined in a spherical trap. The 
aim behind this calculation is to explore whether such a system follows: (a) 
shell-structure, and (b) Hund's rule, in analogy with harmonically trapped 
electrons confined in a quantum dot. We also note that a study of Hund's rule 
for fermionic atoms confined in an optical lattice was carried out recently by 
Karkkainen et a/.|10|. 




The addition energy, i.e, the energy required to add an extra atom, to an 
iV-atom trap is defined as A/x(iV) = n(N + 1) - (i(N), where fi(N)/fi(N + 
1) represents the chemical potential of an N/(N + 1) particle system. The 
chemical potentials, in turn, are defined as fJ,(N) = E(N) — E(N — 1), where 
E(N)/E(N+1) represents the total energy of an N/(N+1) particle system. In 
our calculations, the total energies were calculated using the UHF approach 
for various values of the scattering length and our results for the addition 
energy for an a = 0.01a x spherical trap are presented in Fig. [2J, for the values 
from iV = 1 to iV = 21 . 

For the range of N values studied here, in a noninteracting model the charg- 
ing energy acquires nonzero values A/i(iV) = fkv x , only for iV = 2, 8, and 
20, corresponding to filled-shell configurations. In an interacting model, how- 
ever, Afi(N) should additionally exhibit smaller peaks at N = 5, iV = 14, 
corresponding to the half-filled shells. If the inter-particle repulsion is strong 
enough to split 3s and 3d shells significantly, we will additionally obtain a 
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peak at iV = 18 corresponding to the filled 3d shell, while the peaks corre- 
sponding to the half-filled shells will occur at N = 13, and N = 19, instead 
of N = 14. Moreover, it is of considerable interest to examine whether the 
Hund's rule is also satisfied for open-shell configurations of such spherically 
trapped fermionic atoms, as is the case, e.g., for electrons in quantum dots|TTj. 
From Fig. [2] it is obvious that major peaks are located at N = 2, 8, and 20, 
while the minor ones are at N = 5, and 14, with no peaks at N = 13, 18, or 
19. The heights of the major peaks are in the descending order with increasing 
N, ranging from 1.003Hu x (N = 2) to 0.908^ (N = 20). Additionally, for 
all the open-shell cases, the lowest-energy configurations were consistent with 
the Hund's rule in that, a given shell is first filled with fermions of one (say 
'up') spin-orientation, and upon completion, followed by the fermions of other 
('down') spin orientation. We note that these results are qualitatively similar 
to the results obtained for spherical quantum dots[UJ. Thus, we conclude that 
for the small number of particles considered by us, the shell structure and the 
Hund's rule are also followed by atoms confined in harmonic traps where the 
mutual repulsion is through short-range the contact interaction. 

We have performed a number of UHF calculations on traps of different ge- 
ometries, and scattering lengths, whose results will be published elsewhere. 
However, we would like to briefly state that as the scattering length is in- 
creased, in several cases the ferromagnetic configurations violating the Hund's 
rule become energetically more stable. This implies that for large scattering 
lengths the UHF mean-field approach may not be representative of the true 
state, and inclusion of correlation effects may be necessary. 



6 Conclusions and Future Directions 

In this paper we reported a Fortran 90 implementation of a harmonic oscillator 
basis set based approach towards obtaining the numerical solutions of both 
the restricted, as well as the unrestricted Hartree-Fock equations for spin-| 
fermions confined by a harmonic potential, and interacting via pair-wise delta- 
function potential. The spin-| fermions under consideration could represent a 
two-component fermi gas composed of atoms confined in harmonic traps. We 
performed a number of calculations assuming both attractive, and repulsive, 
inter-particle interactions. As expected, the Hartree-Fock method becomes 
unstable with the increasing scattering length for attractive interactions, while 
no such problem is encountered for the repulsive interactions. Additionally, we 
performed a UHF study of atoms confined in a spherical harmonic trap and 
verified the existence of a shell structure, and that the Hund's rule is followed. 
These results are in good qualitative agreement with similar studies performed 
on harmonically confined electrons in quantum dots, interacting via Coulomb 
interaction. 
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In future, we intend to extend and improve the fermionic aspects of the 
present computer program in several possible ways. As far as problems re- 
lated to fermionic gases in a trap are concerned, we would like to implement 
the Hartree-Fock-Bogoliubov approach to allow us to study such systems in 
the thermodynamic limit, and at finite temperatures. With the aim of study- 
ing the electronic structure of quantum dots, we plan to introduce the option 
of using the Coulomb-repulsion for interparticle interactions, a step which will 
require significant code writing for the two-electron matrix elements. Addition- 
ally, we also aim to introduce the option of studying the dynamics of electrons 
in the presence of an external magnetic field, which will also allow us to study 
fermionic gases in rotating traps. Finally, we plan to implement the option 
of including spin-orbit coupling in our approach, which, at present, is a very 
active area of research. We will report results along these lines in the future, 
as and when they become available. 
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